% k-μ 参数3d图
clc; clear;

% 参数范围
k_values  = linspace(0, 0.5, 30);   % k 范围
mu_values = linspace(2.6, 4, 30);   % mu 范围
z0_values = linspace(-3, 0, 30);    % z0 范围

% ndgrid 可生成对应 (k, mu, z0) 的 3D 网格
[K3, MU3, Z0] = ndgrid(k_values, mu_values, z0_values);

% 预分配分类矩阵 (1=HC, 2=CH, 3=QP, 4=Period)
dynamics_3D = zeros(size(K3));

% 初始化参数
stateX0 = 0.1;
stateY0 = 0;

nTransient = 500;  % 瞬态去除步数
nLE = 500;         % Lyapunov 指数的迭代步数

parfor idx = 1:numel(K3)
    k   = K3(idx);
    mu  = MU3(idx);
    z0i = Z0(idx);  % 记忆阻器初始值

    % 初始化状态
    state = [stateX0, stateY0, z0i];

    % 丢弃瞬态
    for n = 1:nTransient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end

    % 计算 Lyapunov 指数 (LEs)
    LE_vals = LEs(state, mu, k, nLE);
    positive_LE = sum(LE_vals > 0);
    max_LE = max(LE_vals);

    % 判断分类
    if positive_LE >= 2
        dynType = 1;  % Hyperchaos
    elseif positive_LE == 1
        dynType = 2;  % Chaos
    elseif abs(max_LE) < 1e-3
        dynType = 3;  % Quasi-period
    else
        dynType = 4;  % Period
    end

    dynamics_3D(idx) = dynType;
end

% 散点图绘制 3D 分布
figure('Name','Fig2(a)-like 3D','Position',[100 100 800 600]);

% 将 3D 网格与分类结果展开为向量
k_all   = K3(:);
mu_all  = MU3(:);
z0_all  = Z0(:);
dyn_all = dynamics_3D(:);

% 4 种类型对应 4 种颜色
% 1=HC(棕), 2=CH(红), 3=QP(橙), 4=Period(蓝)
color_map = [0.6 0.3 0;  % 1 Hyperchaos
             1   0   0;  % 2 Chaos
             1 0.5  0;  % 3 Quasi-period
             0   0   1]; % 4 Period

% 根据 dyn_all 为每个点分配颜色
colorData = zeros(numel(dyn_all), 3);
for i = 1:numel(dyn_all)
    colorData(i,:) = color_map(dyn_all(i), :);
end

% 使用 scatter3 绘制 3D 散点
scatter3(k_all, mu_all, z0_all, 20, colorData, 'filled');
xlabel('k');
ylabel('\mu');
zlabel('z_0');
title('3D Dynamics Distribution for MCLM');
axis tight; grid on; box on;
view(135, 25);  % 调整视角
